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Abstract 

We study an expansion method for high-dimensional parabolic PDEs which constructs 
accurate approximate solutions by decomposition into solutions to lower-dimensional 
PDEs, and which is particularly effective if there are a low number of dominant principal 
components. The focus of the present article is the derivation of sharp error bounds for 
the constant coefficient case and a first and second order approximation. We give a precise 
characterisation when these bounds hold for (non-smooth) option pricing applications and 
provide numerical results demonstrating that the practically observed convergence speed 
is in agreement with the theoretical predictions. 

Key words: high-dimensional PDEs; asymptotic expansions; anchored ANOVA; financial 
derivative pricing 


1 Introduction 

High-dimensional partial differential equations play an important role in the modelling of 
many real-world phenomena. This is chiefly the case in finance and economics, where one is 
often concerned with probability densities and expectations of economical or financial time 
series modelled by multivariate stochastic processes. These financial applications are the 
main motivation for the analysis in this paper, but we anticipate that similar methods can be 
useful, e.g., in the context of Fokker-Planck equations for high-dimensional chemical systems. 

The computational effort necessary to solve A-dimensional PDEs numerically with stan¬ 
dard grid-based methods grows exponentially with N, a phenomenon referred to as the “curse 
of dimensionality”. Even more sophisticated PDE methods tailored to high-dimensional ap¬ 
proximation, such as those based on sparse grids (see [2] for a survey), are typically not 
able to deal with practical problems where N exceeds about five to eight (see [12] for results 
with sparse finite elements and |2n( 129] for the sparse grid combination technique). A sparse 
wavelet method is proposed in [3Tj, which gives almost dimension-independent convergence 
rates for parabolic equations with non-smooth initial data, as they are typical for derivative 
pricing applications, where the Cauchy data are singular with a singularity located in a low 
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co-dimensional manifold. Numerical experiments for Black-Scholes and stochastic volatility 
models are given in m- 

A radial basis function method for multi-dimensional Black-Scholes PDEs is presented in 
[25] . but only numerical results for one and two dimensions are presented. 

It seems clear that if the approximation is based on tensor product meshes, a compression 
of the data has to take place in the solution process. An active area of research is concerned 
with low rank tensor approximations of such functions (see, e.g., the survey paper [6]), and 
in particular the representation in the so-called tensor train format. The method can be 
targeted specifically towards high-dimensional diffusion problems, motivated by applications 
in finance, and m derive theoretical bounds on the rank of such approximations, which 
break the curse of dimensionality. We are, however, not aware of any successful practical 
application. 

There is a sizeable and growing body of literature which deals with the feasibility of 
integration (cubature) in high dimensions. Several of the ideas from cubature are transferable 
to the solution of PDEs, especially in the situation where the PDE solution can be expressed 
as a high-dimensional integral via a known Green’s function (typically the transition density 
of an underlying stochastic process). 

We will use a PDE-based analysis in this paper. This is to allow for future extensions 
to situations where a Green’s function is not known, or even non-linear and free boundary 
problems, which are not covered by the present analysis. Although the setting in this paper 
- i.e., of linear second order parabolic PDEs with constant coefficients - is in the intersection 
of problems where PDE and cubature methods are applicable, we provide novel error bounds 
even in this case. 

Our approach here is most closely related to cubature based on tensor product decom¬ 
positions. We focus therefore on that segment of the literatureS In particular, [32l [33] 
construct lattice rules by “dimension-by-dimension” integration. This is especially effective 
if there is an ordering of dimensions by importance. Eor instance, |36] investigates the ef¬ 
fective dimension of integration problems in finance using ANOVA decompositions, including 
numerical examples for path-dependent and multi-asset derivatives; [7] combines truncated 
ANOVA decompositions with dimension adaptive sparse grids [5], provides an error analysis 
in mixed Sobolev norms and presents numerical examples for Gollateralized Mortgage Obli¬ 
gations and Asian options. A general framework of decompositions into lower dimensional 
projections is developed in |19j . who show certain minimality and uniqueness properties of 
such decompositions, which contain ANOVA as a special case. 

A similar decomposition is the basis of a PGA-based expansion method for PDEs intro¬ 
duced in [29]. The method exploits the empirical observation that many diffusion processes 
of interest (such as stock prices in equity baskets, or forward interest rates with different 
maturities) have relatively large correlations, which leads to covariance matrices with one or 
a few significant principal components, while the eigenvalues for the remaining eigenvectors 
are an order of magnitude smaller. Transformed into the corresponding basis, the solution of 
the PDE can be reasonably well approximated by a low-dimensional PDE solution restricted 
to the first principal components. Then, dimensionwise corrections (interpretable as approxi¬ 
mate Taylor expansions in the small eigenvalues) can be added to arrive at successively more 
accurate approximations. The key point is that even if the dimension of the original PDE 

^Alternative approaches are Monte Carlo methods or cubature on Wiener space | 23i 1101 I22| . based on 
high-order integration rules derived by exact integration of multivariate polynomials of Brownian integrals. 
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may be very high (say 10 to 100), an approximation with practically sufficient accuracy can 
often be constructed by the solution of a sequence of much lower-dimensional PDEs (say 2- 
or 3-dimensional). 

The relation between PDE expansions and anchored ANOVA is pointed out in [26] and 
[30j : a higher order extension is also sketched in [T3|. A practically useful feature of decom¬ 
position methods is their inherent parallelism, which is exploited in [29| and more recently 
in m- An extension to stochastic volatility models is given in m. where option values 
for baskets with two to eight equities, each with their stochastic volatility process, are com¬ 
puted, albeit with significantly reduced accuracy in the higher-dimensional cases. In [27], 
we demonstrated that these methods can be successfully applied to more complex financial 
market models and derivatives. Eor example, we were able to solve 60-dimensional PDEs for 
Bermudan Swaptions in the LIBOR Market Model with comparable accuracy to standard 
Monte Carlo methods in similar or less (often significantly less) computation time. 

Another line of research concerns the optimal coordinate system as basis for these expan¬ 
sions. PCA-based expansions (PDE-motivated or anchored ANOVA) are well-suited to deal 
with many derivative pricing applications, because market assets typically show high levels 
of correlation, and therefore a transformation to principal components gives a natural order¬ 
ing with rapidly decreasing importance of higher-dimensional contributions. In [14], optimal 
linear transformations of the original coordinates are constructed in order to minimise the 
effective dimension, while [8] proposes non-linear coordinate transformations to extend linear 
PCA. 

In this article, we investigate the theoretical properties of the PCA-based PDE expansion 
method from [29] . We provide a framework for establishing the existence and accuracy of the 
approximate solutions, and give precise theoretical error bounds for first and second order 
versions. While the method has previously been motivated heuristically via Taylor-expansions 
and its relation to anchored ANOVA methods, and its usefulness has been demonstrated by 
successful numerical studies, these are to the authors’ knowledge the first theoretical error 
bounds in terms of the PDE coefficients and the data. 

Specifically, the contribution of this article is to 

• derive error bounds for the PCA-based expansion method for parabolic constant coef¬ 
ficient PDEs, where the initial data have certain mixed-order smoothness; 

• analyse theoretically the applicability to non-smooth data which typically arise in fi¬ 
nancial engineering; 

• demonstrate agreement between theoretical predictions and experimental results. 

The rest of this paper is organised as follows. Section [2] presents the PCA-based expansion 
method for the heat equation and states and discusses the first main theoretical result. Section 
|3| proves sharp error bounds for a first- and second order expansion under sufficient regularity. 
Section 0] provides an alternative construction of the schemes via Taylor expansions and 
discusses the links. Section [5] analyses the applicability to non-smooth data, focussing on 
cases arising in finance. In Section [6l we first generalise the method to constant coefficient 
parabolic PDEs and then give numerical examples which demonstrate agreement between the 
theoretical predictions and numerical results. Einally, Section [7] summarizes the results and 
describes future research directions. 
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2 Set-up and first main result 

To introduce the main concepts, we focus on the A^-dimensional heat equation 


du 

dt 


(1) 

u{z, 0) 

= 9{z), 

(2) 


for z G t G (0,T), A = (Ai,..., Ajv) G M+, and assume Ai > A2 > • • • > 0. 

We will discuss in Section 16.11 how principal component analysis (PCA) can be used to 
transform a general parabolic PDE with constant coefficients to this form, and outline in 
Section [7] how localisation arguments can be used to deal with variable coefficients, with a 
reference to [25] where numerical tests are presented for the variable coefficient setting. 

Definition 1. Given an initial-value problem of the form m and m, and an index set 
j/C{l,...,iV}, define a differential operator 




and an approximation u’^ to u as the solution of 

du 


= O'u'^ 


dt 

= g 


(3) 

(4) 

(5) 


Note that since only operates on the dimensions in the index set u, the problem of 
calculating u’^{zo,T) for a fixed zq G is of spatial dimension |i^|. Calculating the full 
solution u'^{z,T) for all values of z G is still an N^-dimensional problem. 

Definition 2. For a given ^ = {(wi, ni), (w 2 , i' 2 ), ■ ■ ■, (wn, Vn)}, where re* G M and Vi C 
{!,..., N} for all 1 < i < n, we define 

^ w . (6) 

We will refer to as a (truncated) expansion. 

The idea is that ^ encodes an approximate solution for u via expansion into solutions u'^ 
of lower-dimensional PDEs. We will use the notation 

:= — u (7) 


for the expansion error. 

In the following, we introduce the expansions considered in this paper. Consider the solu¬ 
tion u{z, t, A) of ([I]) explicitly as a function of A and define, for some 0 < A*^ = (A?,..., A^) G 
0 < (5A = ((5Ai,..., SXn) € R^, ej the j-th canonical basis vector, 

{Sju){z,t, X^) = u{z,t, X^SXjej), 

{Aju){z,t, X^) = u{z,t, X^ + 6Xjej) — u{z,t, X^), 


4 



such that Sj = I + A,-. Then 


N 


L{z,t,X^ + 6X) = u{z,t,X^) = J|(/ + Aj) u{z,t,X^) = ^ A"u(z,t, A°), ( 8 ) 


ki=i 


ki=i 


ae{ 0 ,l}^ 

where a is a multi-index of O’s and I’s, i.e., the difference operator in each direction appears at 
most with power 1 in each term, and A“ = instance, for N = 2, and omitting 

z and t for brevity, 

u(A? + 5Ai,A^ + 5A2) = u(A?,A^) 


= a;ao« = A(0.0)t 


-|- u{X^ + ( 5 Ai, A2) — 'u(A 5 , A2) + u(A^, A2 + 5A2) — u(A^, A2) 


0 \0 


vO \ 0 ^ 




0 \0 


+ + (5Ai, A2 + (5A2) ~ ^(A^ + (5Ai, A2) ~ 't^(A^, A2 + (5A2) + A2) 




Now consider specifically A*^ = (Ai,..., A^, 0,0,..., 0), 6X = X — A*^, then, for any m > 0, 


N 


“W = EE^““T) + E E^““T) 

j=0 \a\=j j=m+l \a\=j 


(9) 


T (10) 

where Ur^m is seen as an approximation to u, and Ur,m the error, a G { 0 , 1 }^, |a| = 

Then for a G {0,1}-^ one finds that 

A"u(A°) = ^ (-l)"-^u(A° + (5A-/3) if V A: < r : a*, = 0, (11) 

0</3<a 

where is element-wise multiplication, and A“u(A^) = 0 otherwise. The last statement 
follows because SX^ = 0 for all k < r. To demonstrate m, we proceed by induction in |a|. 
The statement is clearly true for a = 0. Consider next a = ej for some j > r, i.e. |q;| = 1. 
Then 


A“u = u(Ao-l-AjCj) — u(Ao) = n(Ao-I-(5A • a) — u(Ao) = (“1)'^ ^u{X^ + 5X ■ (5). 

0</3<a 

Let now 0 < a, 0 7 ^ a with aj = 0 for some j. Then, from the induction hypothesis. 


A“+"^u(A°) = Aj Y, (- 1 )““^u(A° + <5A-/3) 

0</3<o 


= (_i)“-/3n(AO + (5A + A,e,)-/3)- J] (-l)“-^u(A° + 5A •/3) 

0</3<a O<0<a 

Y {-lT~^^^u{X°+ 6X-P)+ Y (-1)““^+^^^(A° + (5A-/3) 

ej</3<a-\-ej 0</5<a 

^ (-l)“+"^-^u(A° + 5A-/3). 

0<p<a-\-ej 

We can make a number of observations: 
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• For A > 0, A“/A“ is a consistent approximation to the mixed derivative in A of order 

a. 

• Ur^m is of the form from Definition [2] for some i.e., a linear combination of solntions 

to (jlHS]). We will use both notations interchangeably. 

• The largest nnmber of different non-zero Xj for any such term is r -|- m, which implies 
that no PDF of dimension higher than r -|- m has to be solved to find Ur^m- 

• Lastly, ([9]) represents an anchored ANOVA decomposition of the function u, for a par¬ 
ticular choice of anchor. See, for instance, [7] or m- 

Example 3. For m = 1, we can write Ur^i, using Definition{l\ as 

N 

Ur,l = _ u{h-,r}^ 

k=r-\-l 

N 

= (l+r-iV)n^^’-’">+ ^ u^F-,r,k}^ (12) 

k=r+l 


such that, with Definition Ur^i = for 




{(1 -Fr-A, {!,..., r}), ,r,r+ !}),(!, {!,,... ,r,r-h 2}), 

(1,{1,,... ,r,r +3}),... ,r, N})}. 


(13) 


We now state the main result on the expansion error for Example [3l We will mainly be 
concerned with smooth solutions from the following class of functions. 

Definition 4. Let 

= {geC^: d{^...di^geC\yi<ii<...<h<N}, 

= {g : M continuous : 3 c > 0 V z G \gi^)\ < c} . 


Then we have the following: 

Theorem 5. Assume g G ( 72 , 2 ,mix ([7}^. Then the expansion error Ur^i satisfies 

d^g 


\u. 




r<i<k<N 


dzldzf 


Proof. See end of Section 13.21 


(14) 

□ 


While published error bounds for ANOVA-type expansions tend to focus on L 2 -type esti¬ 
mates of ANOVA terms (see, e.g., [35] or [38|), we use PDE theory to derive Loo error bounds 
in terms of the mixed smoothness of the components not captured in the expansion. 

Remark 6 . There are no “diagonal” quadratic terms with factors A| and univariate deriva¬ 
tives in This has the important consequence that any solutions which only depend on 

one of the Zk are integrated exactly (as the cross-derivatives vanish), and by superposition 
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any linear combination of such terms. By similar reasoning (and as the second derivative of 
functions is zero), the error Ur,i is zero for initial conditions of the form 

N / N \ 

9i^)= 9k {zi, • • • , Zr, Zf;') + ^ Zi Qi^k {zi, . . . , Z^, Z}f) | , 

k=r+l \ i = r+l / 

i^k 

for any funetions gk and gi^k, i-,k = r + 1,... ,N, i.e., they are functions of zi,, Zr,Zk 
only; hence, the expression (...) is affine in Zj, j ^ {1,... ,r, k}, but arbitrarily nonlinear 
as a function of zi,..., Zr, z^. If an initial condition can be approximated well by functions 
of this form, the expansion error Ur^i will be small, irrespective of smoothness of the data or 
smallness of the eigenvalues and time. 

The smoothness requirement g € (j 2 , 2 ,mix considerably weakened using the fol¬ 

lowing observations. First, for the function from da ED, which only contains the diffusion 
terms relating to an index set z/C {!,...,A^}, we have the Green’s function representation 


u''{z,t) = [ ^''{y,t)g{z -y) dy, 

Jrn 


(15) 


^"(2 /,t) 


n 


exp(-z/i/(4Afct)) 

\/47rAfct 




(16) 


where 5{-) is the Dirac delta function. We can therefore replace the initial condition g in the 
estimate (1141) with a function G which has been smoothed by the heat kernel in the first r 
directions, namely, for u = {1,... ,r}, 


G{z, f) = u^{z, t) 



exp(-?/^/(4Afct)) 


-yi,...,Zr - yr,Zr+i,... ,Zn) dyi... dyr.{17) 


The smoothness of G is all that is relevant for the expansion error. Thus, even if g is only, 
e.g., piecewise smooth, it is only in degenerate cases that G is not smooth everywhere. We 
analyse this in Section [S] and, in particular, give illustrative examples in Appendix [Bj 

Corollary 7. Assume G{-,t) G <72,2,mixed Q (T1\)). Then the expansion error Ur^i 

satisfies 


r<i<k<N 


d^G 

dzldzf 


(18) 


Proof. We first note that for fixed t, the solution v{x,t';t) of 


dv 

v{;0) 


N 


E 

k=r-\-l 


d‘^v 


Gi;t) 


agrees with u for t' = t, v{-,t;t) = u{-,t). Moreover, from (fT^ . 


N 

Ur^l = (1 -|- r — N)v^^ -|- = Uo,l, 

fc=r-|-l 
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where is the solution to 


~w~ ^ 

uW(-,0;t) = G(-,t). 

Application of Theorem [5l with replacements u^v,g^G,N^N — r, Xi^ Xi+r, r —>■ 0 
gives the result. □ 


3 Error bounds for first and second order expansions 

Here, we provide a proof of Theorem [5] and its second order extension, Theorem [TTl Specif¬ 
ically, we derive error bounds in terms of the eigenvalues Aj, time t, and derivatives of the 
initial condition. We first state some auxiliary results. 


3.1 Equations for the expansion error 


We start by formulating PDEs for the expansion error, which will allow us to use standard 
PDE estimates to bound error terms. 


Lemma 8. The expansion error from (EQj is the solution of 

d 


dt 


u 


^ = C'iu^ + f mM^x(0,T), 


#(•, 0 ) = 


w — 1 \ g in 

Jw,u)£^ 


sN 


with source term 


(w,i^)G( 


where C {1,..., A^} arbitrary. 
Proof. See Appendix lA.ll 


(19) 

( 20 ) 


□ 


Although Lemma [8] holds for all index sets we have in mind i' for which 

we will apply it later. We will use the following simple result repeatedly. 

Lemma 9. Let f & and let u be a classical solution to du/dt — Cu = f, rt(-, 0) = 0, where 
C is any second order linear elliptic operator (possibly degenerate), with no zero-order terms, 
and I/I < Ct^ everywhere, for some constant C > 0 and p > 0. Then 

|ri| < — 

p + l 


Proof. The right-hand side above is a super-solution. 


□ 







3.2 Error bounds for a first order expansion 

Consider the first order approximation Ur^i from (|lUp . which has the expansion ^ from Example 
[3l We will now derive bounds via the PDE for the error itself, as introduced in Section [2j 

Lemma 10. For the expansion given in equation dig)) , the expansion error from & is the 
solution of 


with source term 





/ 


Proof. See Appendix IA.2I 


k=l ^ 

0 in , 


E 


fc=r+l 


E 

H 





( 21 ) 

( 22 ) 


(23) 

□ 


We can now give a proof of Theorem [5l 

Proof (of Theoreml^. Set Then, taking the difference between p!]) 

and dH), 


dt 


'y A— _ yx—u 

V x.Eu{p-xm_ y x—u- V A— 

E .E 


u 


{l,...,r,k} _ 


U 




i = r+l 
i^k 


E 




02 

dzf 


N 




E 


i = r+l 
i^k 


52 

^dzp 


(24) 


for all {z,t) E x (0,T], with zero initial condition. Gleaning at (|23l) . we need the second 
fc-derivative of this is itself the solution to the inhomogeneous heat equation 


dt 


( 52 


{l,...,r,k} 


E 



A_ 

dzl 


u 




N 

E 

i = r+l 
ijtk 


54 


' dzldzf 


u, 


(25) 


as is seen by differentiating (IMll twice each with respect to Zi and Zk- In turn, by differentiation 
of dll) twice with respect to Zfc, each of the mixed u derivatives on the right-hand side of (12511 
satisfies a heat equation with zero right-hand side and initial condition 

d^g 

dzldzf ’ 
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so that the norm of the solution is bounded by the norm of this initial condition (this follows 
directly from the Green’s function representation (1151) of the solution to this PDE or the weak 
maximum principle for parabolic PDEs, see e.g. Theorem 8 in Chapter 7 of my Applying 
Lemma [9] twice gives the result. 

□ 


3.3 Error bounds for a second order expansion 

We now consider the second order expansion, ^^,2 from (j^. 

Example 11. The second order expansion has the form Ur ^2 = where 

e = {(l + (iV-r)(iV-r-3)/2,{l,...,r}), 

(2 - (AT _ r), {1,... , r, r + 1}),... , (2 - (iV - r), {1,..., r, iV}), 

(1,{1,... ,r,r+ l,r+ 2}),(1,{1,... ,r,r+ l,r+ 3}),... , 

(1, r,iV-2, iV}),(l, r,iV-l,iV})}. (26) 

Remark 12. Computing according to i26\} requires the computation of one r-dimensional 
PDE, N—r of (r+l)-dimensional PDEs and {N—r){N—r — 1)/2 of {r+2)-dimensional PDEs. 
Compared to the first order scheme, this adds one dimension to the highest dimensional PDEs 
and increases their number by a factor 0{N — r). In return, the approximation order increases 
by one. 

Lemma 13. For the expansion given in equation i26\) . the expansion error ifi from & is a 
solution of but now with source term 


N N N 


f - - ^ 


Q2 q2 , 


. dzfdzl - 

fc=r+l l = r+l i = r + l * 

l^k i^kyl 




(27) 


where vP is the solution to 


diP 


dt 


-C^u’^ = u, u‘'{-,0) = 0. 


(28) 
Ui 

Proof. See Appendix IA.3I □ 

Theorem 14. Assume g G (j‘ 2 P,imx (QHl. Then the expansion error Ur ^2 satisfies 

d^g 


|fir,2(') t) II QQ ^ t 'y ^ PiPjPk 

r<i<j<k<N 


dzfdzjdzl 


(29) 


Proof. Set, for notational brevity, 

(9^ . 


V .— Pi^k 


dzf dzl 


u 




dzf dzl 


u 




such that V '.= w — V \s a, single term in the sum on the right-hand side of (|27p . Then 
differentiating (|28|) with respect to Zi and Zk twice, for = {1 ,... ,r, k} and u = {1,... ,r, k, 1}, 
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respectively, v and w are seen to satisfy 


1 -"“. 

.- ^‘^45”' 


(30) 

(31) 


Taking the difference between p0]l and (f^ gives an equation for v, 

.- 4- 

By differentiation of ([301) twice with respect to zi, it is also seen that 


d 


q2. 




92 , 


q2 q2 q2 

' dz? dz} dz? 


Hence we deduce in turn, using Lemma [9l 


A; 


d'^v 

dzf 




32 q 2 q 2 


< 


< 




dzf dzf dzf 

d2 q2 q2 


(32) 


(33) 


dzf dzf dzf 


N 


N 


X] X] X] 

k=r-\-l l = r-\-l i = r-\-l 
Iz^k i^k^l 


d2 q2 q2 


dzf dzf dzf 


(from (p3]) L 
(from (|32])L 

(from Lemma fT3|) . 


□ 

Remark 15. By extension of Remarkl^ there are no univariate or bivariate quadratic terms 
with factors A| or XfXf in /flP)) . Hence, any solutions which only depend on one or two of 
the Zk are integrated exactly, and by superposition any linear combination of such terms. By 
similar reasoning to before, Mr ,2 is zero for initial conditions of the form 

N N 

g{x') — 'y ^ 9j,k {Xj ) ^k ) 3^1, . . . , Xj.) + 'y ^ Xi gi^j^k{Xj 1 Xk, Xi, , Xr), 

j^k=r-\-l iJ,k=r-\-l 

for any functions gj^k and gij^k- 

Remark 16. One can again weaken the smoothness requirements on g, as discussed in Corol¬ 
lary for the first order case. 

Remark 17. Comparing how i29\) emerges from ^21^ to how m emerges from /fill) , we 
conjecture (but do not prove) for a higher order expansions of the form (0) that 

\\u-Ur,m\\oo < ^ (A - A°)"||^_ 

|q:| = m + l 
aiG{0,l} 
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4 Analysis and construction of approximations by Taylor ex¬ 
pansion in A 


This section provides an alternative derivation of the expansions. This establishes a link 
between the ANOVA view and the literature motivated by Taylor expansions, [291 [I3]. It 
will show that higher order approximation of the Taylor terms is not always advantageous 
(Section [42]) and allows us to show sharpness of our previous results fSection 14.,'ll) . 

Let u{z,t,-) G C"*, i.e., m times differentiable in A, and fix G Using multi-index 
notation, we write D^u = = (^i “ ' • • •' (-^Af “ Ao,Ar)“^. 

Then, by standard multivariate calculus (see, e.g., [18]). 


u{z,t,X) = 

with a remainder term Rm- 


_ f\ _ \0\Q! 

- ^Au(z,t,A°) + Rm(z,t,A,A°,u), 

k=0 |q:|=/c 

If u(z,t, ■) G then an explicit form can be given as 

= (A-AO)“i?“(z,t,A,AAu), 

|Q;|=m-|-l 


I I /* ^ 

R^(z,t,A,A^,u) = -^ / (1 — s)l“l“^(iA?u)(z, t, A°-|-s(A — A°))(is. 
al Jo 

We will analyse the existence of these A-derivatives in detail in Section [5l 


(34) 


(35) 

(36) 


4.1 The first order case (see also Example [3] and Theorem [5]) 

For Ofc := (0,..., 0,1,0,... , 0), where the /c-th entry is 1, a first order finite difference ap¬ 
proximation is given for 5Ak > 0 by 


5Ak 


A"'=u(z,t,A°) = 


i(z, t, A® -|- dAfeCfc) — u{z, t, AO) 


(37) 


where is the /c-th unit vector. Now choose 1 < r < V and set A^ = (Ai,..., A^, 0,... , 0), 
6A = A — AO, and m = 1. Furthermore, assume u{z,t,A) is twice continuously differentiable 
in A. Then 

u{z,t,A) = D^u{z,t,A^) + D'J^u{z,t, A^) + Ri{z,t, A, A^,u} 


|a|=l 


a\ 


N 


= u[z 


, t, A^) -|- ^ Afc- 
fc=r-l-l 
N 


L{z,t,A^ + AfcCfc) - u{z,t,A^) - A|i?^“''(z,t, AO -k Akek,A^,u) 




N 


+ ^ A,A 0 ,m)- k 

k=r-\-l 


Y, AfcAzi?“'=+“'(^,LA,A0,n). 

k,l=r + l 
k^l 


12 








Inserting the finite difference approximations for the first derivatives, 

N 

u{z,t,X) = rir,i + ^ 

k,l=r + l 
k^l 

N 

+ ^ A, A°,u) - A° + Afcefc, A°,n)) . (38) 

k=r+l 

The order (in A) of the second term in (1381) agrees precisely with the results in Section [T21 and 
particularly with Theorem [5j We will analyse in Section [5] the relation between the regularity 
in A and the terms in (jl4[) . The last term is a higher order term for smooth (in A) but 

it is this term which prevents us from directly deducing the compact form of Theorem [5] with 
only first order mixed terms. 


4.2 Alternative difference stencils in A 

The methods above use the standard one-sided finite difference approximation 

Aku{z,t,X^) d . , .Ox , ^ 

-sy-= ) + 0(«A»). 

and, for (5A > 0 and a G {0,1}^ such that 6Xk > 0 if and only if = 1, then 


A'^u(z, t, A*^) „ , n, \ 

--= Dxu{z, t,X ) + ^2 0{akSXk). 

k=l 

Considering (j34p . the question arises if there is any benefit in evaluating the partial deriva¬ 
tives iA" by higher order finite difference approximations such that 

P“n(z, t, A°) = D'^uiz, t, A°) + 0{6X^), (39) 

where is a set of multi-indices containing the orders of all error terms. So for from 
before, Sa = {e^ : > 0}, the A:th canonical vector. For example, [H] proposes to use 

high order compact finite difference stencils introduced in |21j . 

Combining equations (l3^ and (l3ll) gives 

^ _ (\ _ xO'io; 

u{z,t,X) = 

/c=0 \a\=k 
m 

+ E E E ^ + Rrniz, t, A, A°, n). 

k=l |q:|=A; fi&Sa 

If we choose (5A = A — A^ and assume u G in A, the leading order errors are of size 

m 

Y, O ((A - A»)“) + E E E O ((-^ - ^“)“+'’) ■ 

|a|=m+l k=l\a.\=k 
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The order of the first sum is determined by the expansion order (assuming sufficient smooth¬ 
ness) and will be the limiting factor to the overall order in practice. The order of the terms in 
the second sum can be increased by choosing higher order stencils, provided enough smooth¬ 
ness in A, i.e., derivatives up to order a + j3 exist. This comes at the expense of additional 
computational complexity. More specifically, computing 'D°‘u{z,t, X^) requires the values 
u{z, t, A') for different points A', depending on the finite difference scheme. As the dimension¬ 
ality of the PDEs is the same, the complexity increases by a constant multiplicative factor 
and therefore not crucially. All these schemes can be translated into an expansion using the 
z^/^-notation introduced in Section [2l 

Remark 18. A disadvantage of these higher-order steneils is that additional (univariate) 
terms appear in the error and the exaet computation of certain funetions with low 

superposition dimension diseussed in Remark 0 is lost. A similar eomment extends to the 
second-order ease. For the same reason, it does not seem advisable to try and inerease the 
aceuraey of the finite differenee approximations by using a smaller step size than Aj — Ao,i, or 
to eompute the derivatives more or less exactly by, say, algorithmic differentiation. 


4.3 Sharpness of the bounds 

Here, we will show that the error bounds (|14p and (12911 are ‘asymptotically sharp’ in the 
following sense: We will give initial data g such that 

||n-n,,^||oo = (A-A°)“||Il2«^||oo + o(t’"+'|A-A°r+'). (40) 

\a\ =m-\-l 


Take without loss of generality r = 0, A*^ = 0 and 


N 


9{z) = n ^os{zk), 


k=l 


such that Halloo = 5 ( 0 ) = 1 and ||iA^"(7||oo = \D‘^'^g{0)\ = 1 for all a. The solution to ([T] [2]) is 
then 


By (l3ll) one gets that 


N \ N 


i{z, t) = exp -t ^ Afc cos{zk). 


k=l / k=l 


N 


exp 


k=l 


N 


N 


“ 1(e - Ij-FRi 


k\ , 

k=0 \ k=l 


N 


N 




k=l 


—tAi. 




k=l 


i iA/j tXj _ tX]^ _ tA_ 


/ _ e '-'fc - e -\-l] + R 2 


etc, where 


Rm. — 


^ (A-AY^“(z,t,A,A°,n), 


|q:| =m+l 
aiG{0,l} 
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i.e., Rm contains only the mixed terms of Rm in (1351) . One calculates explicitly from (f36]l . 

I I ^ 

f,a ^ m (i_s)l«l-itl«l^(2,t,A0 + s(A-A°))ds = tl“lu(z,t,A°)+ o(tH). 

a! Jo 

From this the claim (jdOh follows. In fact, we see that locally, in a neighbourhood of {z, t) = 
(0,0), (1401) will generally describe the behaviour of the error for smooth g. 

5 Regularity in A and z 

In this section, we establish conditions that guarantee the existence of partial derivatives in A, 
as needed for the Taylor expansions in Section 01 We also derive bounds on the size of those 
derivatives, so that an explicit upper bound for the remainder term in the Taylor expansion 
can be given, thus bridging the analyses from Sections 13.21 and 13.31 on the one hand, and 0] on 
the other. 

In particular, we examine when D'^u exists for non-smooth initial conditions. It is a 
well-known phenomenon that the heat equation has a smoothing effect on non-smooth initial 
data. This has been observed in [9] from an ANOVA decomposition perspective, where the 
authors are, as here, motivated by option pricing problems. Here, we analyze different types 
of payoffs in more detail. Thus, we will be able to answer the regularity question positively 
for a wide class of practically relevant problems. This includes initial conditions which are 
not differentiable everywhere — for example ones with kinks or jumps — as long as the 
non-differentiabilities are localised. Examples 1241 and 1251 illustrate this effect. 

5.1 The smooth case 

Let us first consider the existence and form of partial derivatives D°u for smooth g. 

Proposition 19. Let u be a solution to [l\) with g G (^ 724 .™^ Qjid (j) = then for all 

z € i = 1,..., A^, 

-^u{z,t,X) = tf ^{y,t,X)^{z-y) dy. (41) 

oXi J^N dzf 

If, moreover, g G then for i ^ j 

Q‘2 r Qiq 

(42) 

Proof. This follows directly by differentiating the Green’s function representation (jl5p . □ 

An important point to note from (021) is that as Xi,Xj —>• 0, the integral tends to the 
second mixed derivative of g at z. These are precisely the error terms found in (jl4h . 

5.2 The piecewise smooth case and data smoothing 

The motivation for the research in this paper are applications in derivative pricing, where the 
initial data are typically non-smooth. The explicit examples in Appendix [B] show that this 
does not necessarily mean that the results from Section 15.11 are not applicable for piecewise 
smooth data, in fact, with the exception of degenerate cases, which we will characterise in 
this section, the same convergence orders as in the smooth case hold. 
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Corollary 20 (to Proposition I19p . Let u be a solution to m and G from Assume 

the conditions imposed on g in Proposition [13 on are replaced by the same conditions 
on G{zi, Zr, t) on for fixed zi,..., Then the result still holds for all 

i > j > r and all X> 0 such that Ai ,... ,Xr >0. 

In the following, let X-j = {xi,... ,Xj-i,Xj^i,... ,X]\f) and let dx-j be shorthand for 
dxi... dxj-idxjj^i... dxN- We also write \/A = (\/W) • • •; V^n) and x-y = {xiyi,... ,X]sfyN) 
for the element-wise product. 

Proposition 21. Let g ^ G^ be such that, for 

Gj{z-j,Zj,t,X) = / ^^^’■■■’^~^^x-j,t, l)g{z -VX- {xi,... ,Xj-i,0,Xj+i,.. .,xjv)) dx-j, 

jRN-l 

and every fixed Z-j,t, X, Gj{z-j, ■,t, A) is piecewise C^, i.e., there is a function Gj E (j‘^X,rmx^ 
such that 

N N 

Gj (zj ) = Gj (zj ) ^ Cij max(zj - Uij , 0) ^ bijH{zj - atj ), (43) 

2 = 0 2 = 0 

where H is the Heaviside function and aij,bij,Cij given. 

Then the following holds: 

1. -^u{z,t,X) exists for all z and all X with Xj > 0; 

2. for a given Zj, 


A. 

dX. 


u{z, t, A) 


= t 


X-G, 

dz-j 


A) 


if the derivative on the right-hand side exists. 


Proof. This follows by studying the terms in (|43l) individually. The smooth part is covered 
by Proposition [T^ the piecewise linear Lipschitz and step function parts are as in Examples 
[25] and [24] in Appendix [B] For a more explicit proof see [37]. □ 

Proposition [2T] gives a sufficient but indirect condition for the existence of A-derivatives 
in terms of Gj. The same degree of smoothness is needed for the error bounds in Section 
[3l To understand what conditions on g guarantee that Gj satisfies the requirements, it is 
instructive to consider first the two-dimensional case. The key idea is that integration in one 
direction, zi, helps smooth out discontinuities in Z 2 , as long these are not orthogonal to zi. 

Example 22 (Mixed smoothness). Consider Figure {Jl Out of the five marked points, g is 
twice continuously differentiable only in points p 2 , Ps and p^. Once integrated over zi, the 
resulting function Gi is twice continuously differentiable in p 2 , Pi and p^. Integration helps 
to smooth out the discontinuity in ps. It does not help for pi, where the kink is orthogonal to 
the Z 2 -direction and the discontinuity is orthogonal to zi. Integration over zi decreases the 
smoothness in p^, since g is differentiable in p^ but not differentiable on {(zi,Z 2 ) : 5 < zi < 
10, Z2 = 5}. 
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Figure 1: Schematic of the payoff function for Example 1221 to illustrate how integration in zi 
affects the differentiability in Z 2 ' 5=1 inside the shaded area and 5 = 0 everywhere else. 


Example 23 (Single kink or discontinuity). This is essentially the setting also studied in 
Assume now that 5 : —>■ M is twice continuously differentiable everywhere but on a set 


{( 2 : 1 , ^ 2 ) : f{zi,Z 2 ) = 0 }, 

given by the roots of a smooth function / : M. For simplicity assume a single kink or 

discontinuity. More specifically, we assume that ^ >0 (or <t)) everywhere. 

Then for every value of Z 2 there is at most one value of zi such that f{zi,Z 2 ) = 0, and, 
using the implicit function theorem there exists a smooth function 6 : M —)• M describing the 
location of the kink or discontinuity. 


f{b{z2),Z2) =0 VZ2 G K- 


Writing out the partial derivative now leads to 

d f 

Gi{z,t,X 2 ,X) = ^{xi,t,l)g{z - VX-x)dxi 

dz2 Jr 

•b(z2) roo 

^{xi,t,l)g{z — vX ■ x)dxi + / <h(xi,t, 1 ) 5 ( 2 ; — vA • 

-00 Jb(Z9.) 


A 

dz2 


d 

dz2 

db 

dz2 


{Z 2 )^{h{z 2 ),t,l) g{z - y/x ■ {b{z 2 ) , X 2 )) - g{z - y/X ■ {b{z 2 )'^, X 2 )) 


+ 


[ ^{xi,t,l)-^^(z — y/X ■ x)dxi + [ ^{xi,t,l)-^^{z — ^/X ■ x)dxi. 

J-00 OZ2 Jb(z2) ^^2 


Here, b is a smooth function and 5 is W in Z 2 everywhere but at the kink. Thus every term 
in the sum and the sum itself are continuously differentiable in 22 • This can also be seen by 
writing out the second derivative. 

It is straightforward to generalise this argument to finitely many kink points per Z 2 by 
observing that f is then locally invertible. 
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In higher dimensions, for Zj-derivatives of G to exist in a point it is 

sufficient that either 


g is sufficiently 2 ;j-differentiable for all points in a plane through orthogonal to the 


Zj-axis, or, 


• smoothing of any discontinuities occurs through integration over zi,... ,Zr, which hap¬ 
pens if the location of the discontinuities is not parallel to all zi,... ,Zr. This is guar¬ 
anteed if the set is locally described by 


fizi, ■ ■ 


Zri 1 - 1 ) ' 


.0 .. .0 


) ^n) — 0) 


for hxed ..., z^j_i, Zj_^_l ,..., z%, with smooth /, and 

f (-^1 1 ‘ ‘ ‘ Zri '^r+1) * * * ) 1 ^3 1 ■ ■ ■ ) ) 7“ ^ 


for at least one 1 < k < r, i.e., the kink set is not parallel to the plane and smoothing 
occurs across it in at least one variable z^. 


6 Numerical computations 

6.1 Definition of the model problem 

Consider the Black-Scholes model with stock price processes Si,..., following 

dSi^t = rfSi^t dt + (TiSi^t dWi^t 'i- <i < N, (44) 

where ry G M is the risk-free interest rate and a G are the constant volatilities. For the 
standard iV-dimensional Brownian motion W we take a correlation matrix 


/1 

7 

7 • 

. 7 ^ 

7 

1 

7 • 

. 7 

V 7 

7 

7 • 

• 1 / 


(45) 


for some 7 G (— 1 , 1 ). The covariance matrix S G R^^^ is given by Sjj = aiajpij. In 
the simplest case of identical volatilities ai = ... = aN = a, this implies that there are at 
most two distinct eigenvalues Ai > A2 (with Ai = A2 if and only if 7 = 0) such that A = 
(Ai, A2, A2,..., A2). The normalized eigenvector corresponding to Ai is qi = ( 1 ,..., 1)/i/N 
and the other eigenvectors can be chosen as any orthogonal basis of the orthogonal complement 
of {o’!}. Writing the covariance matrix as a'^{'yE -|- (1 — 7)1), where I is the N x N identity 
matrix and E the matrix full of ones, and noting that Eqi = Nqi and Eq2 = 0, one finds 
the eigenvalues Ai = cr^((A/ — 1)7 -|- 1 ) and A2 = cr^(l — 7). In particular, A2 and the higher 
eigenvalues do not depend on N. 

This covariance matrix is clearly stylised, however, numerical tests with more general 
models confirm that the dominant effect comes from the spectral gap after the rth (here, 
first) eigenvalue, whereas the more detailed structure of the correlations is less relevant. For 
instance, [27] uses a correlation matrix with elements pij = exp(—a|i — j\) for some a > 0 
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in the context of the LIBOR market model, and analyses numerically the impact of 7 on the 
eigenvalues and the numerical accuracy of truncated expansion approximations to Bermudan 
swaptions with N = 5,11,21,41, • • • ,101 tenors. Closer to the test cases here, [2U] gives 
results for arithmetic basket options, but with an empirically estimated covariance matrix 
for the DAX, i.e. N = 30. There, the eigenvalues jump by roughly a factor of 10 after the 
first, and then decay slowly but exponentially. A similar size jump is observed in Figure [2] for 
7 ~ 0.4 to 0.5, a typical correlation between equity prices. The accuracy for the example in 
[29] is 0.73 basis points. 

For simplicity, we choose rj = 0. The arbitrage-free price V{Sq, 0) at time 0, of a financial 
derivative with payoff hlSj^) at time T is then given by 

ViSo,0)=E[h{ST)\So]. (46) 

We can estimate this value with desired accuracy using, for example, Monte Carlo simulation 
with a sufficiently large number of samples. This provides a reference to compare the PDF 
solution to. 

The expectation (l46]) is given by the solution u(log(S'o), T) of 

du _ d'^u ^ aj du 

dt ^ dxidxj ^ 2 dxi' 

u{x,0) = ff(exp(x)), 

where we have used a similar (logarithmic) transformation as in Appendix iBl 

For constant, positive semidefinite coefficient matrix S = {(7iajPij)i<ij<N, we can choose 
Q to be the matrix of eigenvectors of S sorted by eigenvalue size, i.e., 

Q = {qi,...,qN), ^T.qi = Xiqi, Ai > ... > Aat > 0 , (47) 

and get A = diag(Ai,..., Aat) a constant diagonal matrix. 

Introducing the new coordinates 

z{x,t) = {x + fit), 

with fi = and inverse transform x{z,t) = Qz — fit, leads to the heat equation 

(Hill]) with g{z) = h{x{z,0)). 

6.2 Numerical approximation of low-order terms 

For each low-dimensional term in the decomposition, a PDF on an unbounded spatial do¬ 
main is solvedH To avoid the introduction of artificial boundary conditions necessary when 
localising the domain, for each coordinate Zi we map the interval (— 00 , 00 ) to ( 0 , 1 ) via 

1 /, ^ 1 
yi = - arctan {piZi + Ci) + -, 
vr z 

with parameters hi and c*. 

^In the present case of constant coefficients, semi-analytic or Fourier methods would be even more efficient 
to solve these sub-problems. 
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Under a standard growth condition on the solution at infinity, the resnlting PDE is fully 
specified without boundary conditions at Zi G {0,1}, because the resulting non-constant 
coefficients of the transformed diffusion equation vanish sufficiently fast at the boundaries 
(see [29l|39]). For options with unbounded payoffs towards infinity, such as European calls, 
we introduce an artificial numerical cutoff. This does not impact the computed option value 
noticeably. 

We consider an equidistant mesh in the yi coordinates, such that in original coordinates 
the mesh is denser in the most relevant regions. This density can be adjusted by choosing 
the parameters bi and Cj, which is performed heuristically. 

We use the Crank-Nicolson time discretization scheme with central spatial differences. In 
the one-dimensional case, solving the PDE is then straightforward. In the two- and three- 
dimensional case, we combine this with an Alternating Direction Implicit (ADI) factorisation 
as detailed below, such that the resulting tridiagonal matrix systems can be solved efficiently in 
linear time (i.e., proportional to the system size). An initial LU factorisation of the tridiagonal 
matrices gives significant further speed-up. 


6.3 Variance reduction of Monte Carlo estimators 

We consider now a simulation-based method for solving the PDE ([ll [2]), making use of the 
expansion method as outlined in Section [2l We will use this for numerical verification only, 
but envisage this to be a good standalone variance reduction method more generally. 

We recall from ([8]) the decomposition 

tt(2,0,A) = 5]] 5A“A“u(z,0,A°). (48) 

We consider independent estimators Ua for in (|48]l . and an estimator 

U = ^ Ua 

oe{o,i}^ 


for u. If the single sample variance of Ua is fa, then, using Na independent samples for Ua, 


Var 



= E 

ag{0,l}^ 


Vg 

Na' 


using independent samples for different a. We expect from the theoretical results in Section 
[3l given enough smoothness, Va = O(A^) for the first order corrections, and Va = O(A^A^) for 
the second order corrections. 

To achieve an overall mean-square error of size e for the first order terms, one thus needs 
0(A|e“^) MC samples, instead of the 0{e~^) required for the naive approach of simulating 
U directly with a single set of samples. In the cases below, Xk was often in the order of 10“^, 
which resulted in a significant run time reduction by a factor of about 10, 000. 

Estimating the first order approximation ui^i in that way, requires N individual estima¬ 
tors. Computing the whole expansion in this way would require 2^ estimators and therefore 
would be infeasible for large N, and would stop to be advantageous already for moderate N. 
A better approach is to compute an accurate approximation to ui^i using finite difference 
methods, and then use a Monte Carlo estimator U\^\ for uip, the difference between u and 
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Figure 2: Absolute difference A (marked x) between PDF and MC results with 3a error 
bounds (-), for varying A 2 and payoff weights ui, 002 and wa from (jM]), (fSOl) and (ISTI) . using 
the first order expansion method. The best fit exponents for the differences are 2.04 ± 0.11, 
2.03±0.11 and 2.21±0.17 (95% confidence bounds from a linear regression on a log-log scale). 


ui^i, with a much smaller variance than estimating u directly, and similar for higher order. 
This is the approach taken in the computations in Sections 16.41 and 16.51 for the smallest values 
of A. 


6.4 Convergence of the first order expansion 

We consider hrst a European arithmetic basket option with maturity T = 1 and payoff 


9{S) 


N 


max 


'^LOiSi - K,0 


\i=l 


where u G is a vector of weights and K is the strike, which we set to K = 100. We study 
the solution at the coordinates Si^ = 100 for all i. The numerical tests are for N = 10 and 
ai = 0.2 for z = 1,... , A. Other details are as specified in Section [ 6 Tl 

Because of the structure of S with one dominant eigenvalue, we choose r = 1 and expand 
around the point A*^ = (Ai,0,... ,0). Figure [2] shows results for 

wi = (1/10,1/10,1/10,1/10,1/10,1/10,1/10,1/10,1/10,1/10), (49) 

u,2 = (4/30,4/30,4/30,4/30,4/30,2/30,2/30,2/30,2/30,2/30), (50) 

(ua = (1/4,1/4,1/4,1/4,1/4,1/4,1/4,-1/4,-1/4,-1/4). (51) 


Plotted is the difference A between the first order expansion solution from Example [3l 
and a Monte Carlo estimate of the true solution, such that A can be considered an estimate 
of the error ui^i from (jlOp . with r = m = 1, and Theorem [5j By varying 7 in (1451) . we vary 
A 2 = As = ... = Aiv, and can therefore examine the expansion error as a function of the 
eigenvalues. 

The one- and two-dimensional PDFs for and for k = 2,..., A were solved 

numerically as described in Section 16.21 Specifically, we use the Peaceman-Ratchford ADI 
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scheme [23] with J = 800 grid points in each spatial direction and M = 50 time steps. The 
MC estimator was constructed by exact sampling of the joint lognormal distribution of the 
asset price vector St from (144 p and computation of the sample mean of the payoff ^(S'r) over 
10^° samples to find 14(5'o,0) by (j36|). 

For very small A 2 as we need for the study of the asymptotic error, the variance of 
straightforward MC was too large to assess the accuracy of the expansion approximation. 
To verify that u — ui^i = 0(A|) by a standard MC estimator for u and PDF solution for 
ui^i, one needs a standard error of ©(A^) or smaller for u, i.e. ©(A^) samples. For the four 
smallest values of A 2 in the tests, A 2 < 0.004 and A| « 2.5 x 10“^°, such that the results for 
a feasible number of samples were too noisy for a meaningful comparison. We thus chose a 
control variate approach as outlined in Section 16.31 in which the difference between the full 
and approximated solution was calculated via a MC estimator using the same paths for both 
terms. 

Figure [2] suggests that the difference A between the truncated and full solution follows a 
power law of the form A ~ A^. The results reported in Figure [2] are consistent with the theo¬ 
retical order of 2 established in Theorem[5l In absolute terms, the error for practically relevant 
correlations around 0.5 is approximately 10“^ (with the exception of the more challenging 
payoff given by ui^), so approximately one basis point in relation to Si^o = 100. 


6.5 Convergence of the second order expansion 


We consider the same arithmetic basket option and model as in Section 16.41 but with N = 5 
(as opposed to A = 10 in the first-order case). This allows us to examine the expansion 
order with significantly reduced computation time, as more (now O(A^) instead of 0{N)) 
and higher-dimensional (now 3-dimensional) PDFs need to be solved in the second-order 
expansion ui ^2 from Fxamplellll 

The one-, two-, and three-dimensional PDFs for, respectively, and 

for 2 < k < I < N where solved numerically as described in Section 16.21 In addition to 
the one-and two-dimensional PDFs already discussed in Section 16.41 we used Brian’s ADI 
method mm with J = 500 grid points in each spatial direction and M = 50 time steps for 
the numerical solution of the three-dimensional PDF. The MC estimator was constructed in 
exactly the same way as in Section 16.41 For the four smallest values of A 2 we again directly 
compute the difference between full and approximate solution via MC, as described in Section 

[Ql 

The matrix of eigenvectors is 


Q 


0.4472 

"k 

kr 

kr 

kr ^ 

0.4472 

kr 

kr 

kr 

kr 

0.4472 

kr 

kr 

kr 

kr 

0.4472 

kr 

kr 

kr 

kr 

0.4472 

kr 

kr 

kr 

^ / 


where the ★ entries can be any set of normalized, orthogonal vectors which span the N — 1 
dimensional subspace orthogonal to the first eigenvector (since A 2 = ... = \n)- For weight 
vectors ui closely aligned with the first eigenvector, such as a; = (1/5,1/5,1/5,1/5,1/5), the 
second order scheme is extremely accurate and the resulting absolute error is considerably 
less than 10~® even for the larger of the values of A 2 considered here. We instead choose 
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Figure 3: Absolute difference A (marked x) between PDF and MC results with 3a error 
bounds (-) for varying A 2 and payoff weights ui and 002 from (IMD and (I^Ut) . using the second 
order expansion method. The best fit exponents for the absolute differences are 3.04 ± 0.28 
and 2.76 ± 0.45. 


payoff weights 


cui = (1,-1,1,-1,1), 

a ;2 = (3/2,3/2,-1/2,-1/2,-1), 


for which the error is relatively large, since they are not aligned with any of the eigenvectors!! 
This allows us to show convergence across a wide range of 7 values (see (j45p ) in Figure El 
The data points again appear to follow a power law A ~ A^. The best fit exponents p from 
Figure El compare well to the theoretical order of 3. 


6.6 Computational times 

In this section we briefly discuss the computational times for the expansion solutions in 
relation to Monte Carlo estimation. Table [H gives CPU times for the set-up in Sections 16.41 
and 16.51 

For the first order expansion (r = 1, s = 1) of Section [6.41 CPU times are dominated by 
the numerical solution of the nine two-dimensional PDFs involved, whereas for the second 
order expansion (r = 1, s = 2) of Section [6.51 the dominant contribution comes from the six 
three-dimensional PDFs. 

The time for the first order expansion with ten assets is roughly one minute, compared to 
five hours for Monte Carlo, and for the second order expansion with five assets, roughly ten 
hours, compared to two hours with Monte Carlo. A few comments are in order. 

3The angle between qi and cai and that between and U 2 are both just under 80 deg, i.e., closer to 
orthogonal than parallel. 
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Method 

MC {N = 10) 

ADI {N = 10, s = 1) 

MC (A = 5) 

ADI (A = 5, s = 2) 

Paths / Points 


M = 50, J = 800 

10 ^*^ 

M = 50, J = 500 

CPU time (sec) 

17902 

68 

7737 

36445 


Table 1: Computational times on a 2.8 GHz Intel Core 2 Duo processor with 8 GB 1067 MHz 
DDRS, running Matlab R2016b. The ADI schemes use M time steps and J mesh intervals 
in each direction, i.e. 6.4 x 10^ unknowns in the 2-dimensional case {s = 1) and 1.25 x 10® 
unknowns in the three-dimensional case (s = 2). 


The ADI computing times are broadly similar to the ones reported in m for two- 
dimensional PDEs and in [H] for three-dimensional PDEs, extrapolating to our finer meshes 
and taking into account minor differences in the set-up. 

The numbers of Monte Carlo paths, ADI time steps and mesh points were chosen to make 
the numerical error much smaller than the expansion error, as the latter is the one of interest 
for the purposes of this paper. In practical applications, there is no benefit in reducing the 
Monte Carlo and finite difference error far below the expansion error, so that a significantly 
smaller number of paths and points can be used, reducing the computing time. 

Specifically, the finite difference error for the above setup is substantially below 10“®. 
For a practically acceptable accuracy of one basis-point, i.e. 1/100 % of S'o = 100, equalling 
10“^, it would suffice to take a quarter of the steps, J = 200 and M = 12, reducing the 
computational time for the first order expansion by a factor of 2® to about 1 sec. Equally, a 
95% Monte Carlo confidence interval with width 2 • 10“^ could already be obtained with lO'^ 
samples, in 18 secs. 

The Monte Carlo estimator samples the joint log-normal distribution in a single time step, 
which would not be possible if a more general model was chosen, which biases the computing 
times by a factor of 12 or 50, respectively, in favour of the Monte Carlo method. So in a 
realistic set-up, given the estimated run-times from the previous paragraph, the combined 
first order expansion finite difference method will be about a factor of 12 x 18 ~ 200 faster 
than a Monte Carlo estimator with similar accuracy. 


6.7 The effect of kinks 

To demonstrate how the solution behaves around kink points, we switch to a geometric basket 
option with = 10 and payoff 


9{S) 


max 




max 



N 


oji log Si 


\i=l 



See Example 1251 for the two-dimensional case. Due to the log-normality, this results in a linear 
kink, for which we can easily choose vectors of weights u which are either exactly parallel, 
orthogonal, or otherwise oriented in relation to the eigenvectors of S to best examine this 
effect. 

For Lj which are not orthogonal to qi, such as uj = qi ot uj = qi + q 2 , this expansion scheme 
is extremely accurate, even at the kink point K = 1. This was expected from the previous 
section. We thus focus on the orthogonal case where no smoothing occurs (see Example [25] 
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Figure 4: Left: Absolute error A (marked x) between PDF and MC solution for a geometric 
basket with payoff weight ivi (bottom) and UJ 2 (top) around the kink point at AT = 1 for 7 = 0.6 
for fixed spot price and varying strike. The option value is around 5 • 10~^. Centre/Right: 
Absolute errors A (marked x) with 3a error bounds (-) for oji with different values for A 2 
at K = 0.5 (centre) and K = I (right). The best fit exponents for the errors are 2.01 ± 0.13 
and 0.61 ± 0.10 (95% confidence bounds from a linear regression on a log-log scale). 


in Appendix |B|), namely0 


= (-0.1160,0.0929, -0.6527, -0.1121,0.6986,0.2091, -0.0438, -0.0758,0.0000,0.000), 

W2 = (0.1130, -0.0607, -0.1708, -0.2057,0.8971, -0.2467, -0.1831, -0.1085, -0.0345, -0.0001). 


The kink point is now located at strike price K = 1. Figured] shows that there is indeed a 
marked increase in the error size around the kink and that the scheme is very accurate away 
from it. The estimated convergence order away from the kink (see Figured]) is consistent with 
the theoretically predicted order of 2 . At the kink, the order of convergence appears reduced 
to below 1, although we expect it to be 1 from Example 1251 

It is worth emphasising that a high order of convergence is achieved in all but the at- 
the-money, orthogonal case, despite the non-smooth initial data. This is in line with, and 
empirically validates, our theoretical analysis in the previous sections. 

7 Conclusion and outlook 

In this paper, we provide a thorough theoretical analysis of the PCA-based expansion ap¬ 
proach for parabolic PDEs first introduced in [29]. Building on our previous work in m , we 
give a description of a method which is general enough to cover a wide range of important 
cases and yet allows concrete error bounds in special cases, such as the ones given in Section 
[3l We illustrate our theoretical results with numerical experiments, which show accuracy and 
convergence in line with theoretical predictions. 

This work focuses on the case of constant coefficients, where the PDE can be transformed 
to the iV-dimensional heat equation via rotation into the eigensystem of the covariance matrix. 
Using PDEs for the expansion error itself, a careful analysis of the applicability of the method 
and the size of the expansion error was possible. 

For non-constant coefficients, one can approximate the PDE via localization. Consider an 

^These two seemingly arbitrary eigenvectors are taken from a MATLAB decomposition for E. 
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A/^-dimensional linear parabolic PDE 


dv 


N r.2 J-y r. 

. . O V ^ , , ov 

» 'J = 1 7 = 1 


(52) 


where T^j = aiGjpij is a function of x and t. For simplicity, we have omitted a zero order 
term, but a term of the form c{t)v can directly be accounted for by the definition of a new 
function exp(Jg c{T)dT)v. 

A transformation to the heat equation like in the constant coefficient case of Section 16.11 
is generally no longer possible. 

A simple approximation is to replace S(x,t) and p{x,t) by S(xo,0) and /u(xo,0), respec¬ 
tively, for some fixed xq. If one is interested only in the solution at (xo,t) for not too large 
t, this may be justihed. One can then proceed as in Section 16.11 This approach leads to a 
localization error which depends on the rate of change of the PDE coefficients in the vicinity 
of (xo,0). 

A concrete example where this approach can be successful is the LIBOR Market Model, 
where the LIBOR rates process has a strongly non-linear drift in the pricing measure, resulting 
in non-constant coefficients for the first derivatives in the pricing PDE. In |27] . the drift is 
“frozen” at the initial value of the process, leading to a constant coefficient PDE, which is 
then approximated by an expansion. The numerical tests for Bermudan and path-dependent 
derivatives illustrate the error from the constant coefficient approximation, as well as the 
expansion error, and indicate that for moderate maturities both are within an acceptable 
range. 

Alternatively, one can adapt the PDE expansion method to directly incorporate non¬ 
constant PDE coefficients. By choosing Q as the eigenvectors of S(xo,0) and transforming 
the PDE ([5^ to the coordinate system z = z{x,t) = Q^{x -|- p{xo, s) ds), one obtains a 
PDE 


du 


1 

2 


N 




d'^u 

dzidzj 


N 

+ '^Ki{z,t) 

i=l 


du 

dzi’ 


(53) 


with coefficients A and k such that Ajj(zo,0) = 0 for z / j and Ki{zo,0) = 0, where zq = 
z{xo,0). 

The dimensionality of the PDE can be reduced by setting the diffusion and drift coefficients 
to zero for indices i ^ v-, for a given index set v. More specifically, we set Ajj = Ajj = 0 and 
ki = 0 for i ^ ly and for all 1 < J < A^, and Ajj = Ajj and ki = Ki otherwise. One can then 
define an approximation u'^(z,t) as solution to (|53p with A and k replaced by A and k. 

The dimensionality of the PDE is then |i^|. The index sets of interest are of the form 
Ur = {!) • • • 5 f}, U {k} for some r < k < N, or UrU {k, 1} for some r<k<l<N, in the 
spirit of this paper. 

It is shown in [28] how to build dimensionwise decompositions in this case and numerical 
results are presented for variable (in time and space) volatilities and correlations. Those 
preliminary tests in [28| show that good accuracy is often already obtained for the first order 
version even in the variable coefficient setting, but for fast varying coefficients the error can 
be significantly larger than for constant coefficients. A sketch of a possible analysis in this 
case is given in [28], but theoretical error bounds for variable coefficients remain an open 
question. 
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It would be interesting to conduct further tests on practically used local or stochastic 
volatility models. For options on a baskets of N assets, where each asset price process 
has its own stochastic volatility, the above method applied directly to the 2A^-dimensional 
pricing PDF can be expected to work well when the asset prices are strongly correlated, and 
strongly negatively correlated with their stochastic volatilities, combined with slowly varying 
volatility. This regime should be close to the constant coefficient case and the eigenvalues of 
the “frozen” covariance matrix will decay rapidly. Interesting adaptations of the method could 
take advantage of the fact that the payoff is only a function of the assets and not the volatility, 
that the inclusion of each stochastic volatility can be seen as adding a correction, and one 
might even combine the eigenvalue expansion with asymptotic expansions with respect to fast 
and slow volatility scales, see [TH] . 

It is conceivable that the method might also be successfully applied to partial integro- 
differential equations arising from multivariate jump-diffusion processes, if a suitable expan¬ 
sion around a dominant component can be defined, e.g., if the jump process can be interpreted 
as a time-changed Brownian motion. 
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A Proofs of lemmas 

A.l Proof of Lemma [8] 

By definition, 

wC'^u'^ — Cu 

dt ^ 

= [£"«-£"«] H wC''u'' + [C''i -C''i]u-Cu 

= jy'i[u^-u]+ w[jy'- jy'i]u''+ [0'^ -L]u. 


A.2 Proof of Lemma 1101 

We use Lemma [51 The relevant differential operators are 


N 


52 


^ = . 


£{!■ ■■■’’'} = y] Afe 


k=l 


dzV 


£{l.....r-.fc} =£{1.-.'-} + A, 


dzV 


and so on. Using this we can rewrite equation (HU) as 





N 

E 

k—r+l 


£{l,...,r,fc} _ £{l,...,r} £{1.£ 


E^'= 



N 

+ Afe 

k—r-\-l 




{l,...,r,fc} 



A.3 Proof of Lemma 1131 

We start with equation dun 


and insert explicitly ^ from (155)) to obtain 


’'■>w« + (l + (iV-r)(iV-r-3)/2) [£{!’■■■’’'>-£{1’■■■’’■> 
ot L 
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+ (2-(fV-r)) Y 


k—r+1 


U 




N N 






k—r-\-l l—k-\-l 


U 






N 


= E^fc 7 j; 2 ^^ + o-(^-’’- 2 ) y]] Afc 


fc=l 

A A 


k—r-^1 


dzi^ 


92 , 92 


TV 


92 


E E + .^ ■'‘M" 

fe=r+l Z=fc+1 ^ ^ ^ ^ fe=r+l ^ 


(54) 
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Observing that the double sum contains exactly {N — r){N — r — 1) terms we can rewrite this as 


= ^Afe^u«-(iV-r-2) ^ Xk^u-{N-r-2) ^ -- 


fe=i 


' dzi 


k=r+l 


-dzl 


k—r-\-l 


N N 


k—r-\-l l—k-\-l 


52 , 92 
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92 
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92 




= Z - 2) E . 


fe=r+l 
N N 


dzl 


fe=r+l 


'5z2 


' dzl 


-dzl 


E E (^4 "^‘ 5 )'“''. 


^ k—r-\-l ^ fc^r+l Z=fc+1 

In other words, solves a non-homogenous r-dimensional heat equation with source term 


N N 


k—r-\-l l—k-\-l 


92 . 92 
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52 




fe=.+i 




- n . 


We know from in the proof of the first order bound that — u is the solution 

to 


•'{l,...,r,fe} _ 

5t 


02 ^ 02 


‘5z2 


dzl 


(55) 


2 e{l,...,r,fc} * 

Using the Green’s function associated with 

pT I- w 

V Xi-^u{z-y,t-s)dyds. 

Jo Js.^ -• 




Similarly, it follows 

dz^ 

i—r+lA^k,l 2 

This allows us to expand the expression for /. We leave out the coordinates to increase readability. 
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{l,...,r,k,l} 


fT r f)2 

(z,t) = - / V \j^u{z-y,t-s)dyds. 

Jo J^^ TT'^i., dz^ 
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94 


9 ^ 25^2 


Both lines contain {N — r){N — r — 1){N — r — 2) individual terms. This allows us to combine them 
into 
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- Y E / / AfcA,[$li. $ 11 . ’’■.'=1] 
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Now from (l28l) 




{z,t) = [ f 

Jo Js.^ 


94 


dzldzf 


= r [ ^^^-^^’^Hy,s)^^u{z-y,t-s)dyds 

Jo jrn ozf^az^ 


94 


and similar for Inserting this in (I5S1) . the result follows. 


B Analytical examples with non-smooth data 


For illustration, we study the two-dimensional Black-Scholes model, where the value function of a 
European option on a pair of assets S' 1 , 5 ' 2 , with payoff h{Si,S 2 ) at time T is given by exp(—r(T — 
t))V{Si,S 2 ,t), where 


0, (56) 

h{Si,S2), 


and CT 12 = p(JiU 2 - 

We exploit the log-normality of the Black-Scholes model by using logarithmic coordinates {xi,X 2 ) = 
(log5'i,logS'2). We now make a specific choice of parameters which simplifies the equations below, 
namely ctJ = <t| = 2r = for some ct > 0. For a general choice of parameters, a further transforma¬ 
tion of (I56|) to eliminate the drift term is needed, but we avoid this complication here as it does not 
change the essence of the problem. Elementary calculation shows that the rotation 



leads to a particularly simple form, the heat equation 

du d'^u 

u(-,0) = g 

for some g, in backward time t ^ T — t. Here, 

Ai=CT^(l-fp), A2=cr^(l-p), 


so that for p close to 1 or —1, i.e., near perfect correlation or anti-correlation, the problem becomes 
close to one-dimensional and the asymptotic expansion can be expected to work particularly well. 

The purpose of the following examples is to show that the optimal convergence order is only lost 
in degenerate cases, which we will describe precisely. Because of the symmetry, we only expand in A 2 
without loss of generality. 

Example 24 (Digital call on geometric basket). Consider the payoff function 


hisi,S2) = l[o,oo)(srs2'- 2) 

with real parameters pi and p 2 , which leads to 


g{zi,Z2) 


^ [0,oo) 



Pi + P2 


-Zl 


Pi - P2 


-Z2 



Note that the level curves of the argument of 1 are straight lines. We consider this as a “local” 
(i.e., linearised) model for the more general case where discontinuities occur along a smooth curve. 
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We will consider three different choices of ni and fi 2 , o-nd analyse the expansion in A 2 , for fixed 


Ai. 

Case 1: = 1, /i 2 = 1 (i.e., gradient of payoff orthogonal to expansion direction). 

In this case, the value of the boundary function depends only on zi and correspondingly 

d 

-^u{z,t,X) = 0 . 

Case 2: yLi = 1, p ,2 = —1 (i.e., gradient of payoff parallel to expansion direction). 

The value of the boundary function depends only on Z 2 . We find that the derivative in A 2 exists and, 
by a lengthy calculation, is given by 

9 22 -log 2 exp(-(z 2 - log2)V4tA2) 

which goes to 0 for A 2 —>■ 0, regardless of the value of Z 2 . In this specific case, the existence of the 
derivative at the kink is due to the symmetry of the PDE and initial condition. 

Case 3: /ii = 2, /i 2 = 0 (i.e., gradient of payoff at 45 deg angle to expansion direction). 

Here, again by elementary but lengthy calculation. 




log 2 - Zi - Z 2 / (log2 - Zi - Z 2 )^A 2 \ 

+ A2^ 4tAi(Ai+A2) / 


In particular, the derivative vanishes at the kink where zi + Z 2 = log 2. Everywhere else, the derivative 
approaches a non-zero bounded value for A 2 —>■ 0 for Ai > 0 fixed. 


Example 25 (Standard call on geometric basket). Consider the same setting as in the previous 
Example \Zf\ but with payoff function 


h{si,S 2 ) = max(s('^S 2 ^ — 2,0), i.e. 

g{zi,Z 2 ) = max I exp 


piE PL2 . Ti- 

—4- 


- 2,0 . 


Case 1: Eor the first case, fii = 1, p ,2 = 1, we have again 

= 0 . 


Case 2: Eor the second case, p-i = 1, /i 2 = —1, we find by straightforward calculus 




^gZ 2 +A 2 i j g—( 22 —log iC+2A2i)^/4iA2 

2 1 V^TrtA^ 


<i) 


Z 2 - logitr + 2X2C 
\/2tX2 , 


(57) 


where $ is the standard normal cdf This derivative exists for all Z 2 G R and all X 2 > 0. For A 2 —>■ 0 
the limit exists and is well-defined provided Z 2 log 2. 
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If Z 2 = log 2, the right-hand side of (153 does not converge for A 2 —!> 0. We can, however, Taylor 
expand in instead of X 2 , since 


d 


d 

i(z 2 , t, X 2 ) = 2y/^——u{Z2, t, X 2 ) 
OA2 


exists in the limit X 2 —^ 0. This expansion will lead to an expansion error of size 0{X2) instead of 
0(A2). Note that the PDE approximation method itself does not have to change. 


Case 3: The third case is analytically lengthier than Example 1^41 o^nd we have therefore omitted it 
here. We refer to Section fST^ and particularly Examples\^ and\^for a qualitative discussion. The 
conclusion will be the same as in Case 3 of Example\24\ i.e., bounded derivative with respect to X 2 for 
X 2 ^ 0 for Ai > 0 fixed. 
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